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We developed analytical and numerical methods to study a transport of non-interacting particles 
in large networks consisting of M d-dimensional containers Ci , . . . , Cm with radii Ri linked together 
by tubes of length Uj and radii Uij where i, j = 1,2,..., M. Tubes may join directly with each other 
forming junctions. It is possible that some links are absent. Instead of solving the diffusion equation 
for the full problem we formulated an approach that is computationally more efficient. We derived 
a set of rate equations that govern the time dependence of the number of particles in each container 
Ni(t),N2(t), . . . , iVjvf(t). In such a way the complicated transport problem is reduced to a set of 
M first order integro-differential equations in time, which can be solved efficiently by the algorithm 
presented here. The workings of the method have been demonstrated on a couple of examples: 
networks involving three, four and seven containers, and one network with a three-point junction. 
Already simple networks with relatively few containers exhibit interesting transport behavior. For 
example, we showed that it is possible to adjust the geometry of the networks so that the particle 
concentration varies in time in a wave-like manner. Such behavior deviates from simple exponential 
growth and decay occurring in the two container system. 



I. INTRODUCTION 

The goal of this work is to find a method that de- 
scribes diffusive transport of particles on a network built 
up of spherical containers connected by tubes. It is pos- 
sible that not all containers are connected to each other 
and tubes may join together forming junctions (without a 
container being present) . In such a way one can generate 
an enormous number of network topologies. One example 
is show in Fig. ^ The networks consists of M containers 
(reservoirs) labeled C\ , . . . , Cm of radii Ri connected by 
tubes of length l^j and radii a.y where i,j = 1,2,..., M. 

Our work is motivated by experiments discussed in 
refs. 0, H and 0, U E|. The first set of references 
describes how to create and manipulate microscale sized 
compartments (vesicles) connected by nanotubes. These 
structures can be applied in a number of ways 0]. For 
example, the non-diffusive (forced) transport was studied 
in P]. In this work we study different kind of transport 
were only passive diffusion is allowed. Also, by including 
the reactions in the theoretical setup one could describe 
biochemical reactions in a milieu close to their natural 
habitat. This idea was pursued experimentally in ref. 



The second set of references deals with networks of 
chemical reactors of macroscopic sizes connected to each 
other by tubes. The exchange of reactants is mediated 
through the tubes and controlled by pumps. It was shown 
that it is possible to use this device to carry out pattern 
recognition tasks. In making the device smaller the ex- 
ternal pumping can be removed and the transport can 
be limited to pure diffusion. This kind of setup is close 
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to the situation studied here. 

For large networks, links connecting opposite sides of 
the network may be rather long. Accordingly, one can not 
expect an exponential decay of the number of particles 
in the containers and this is the situation we are mostly 
interested in. Obviously, to describe such a situation one 
can attempt to solve the diffusion equation numerically 
and obtain the distribution function p(f, t) that describes 
how particles spread throughout the network. 

Obtaining the solution of the diffusion equation for a 
complicated geometry for a large networks gets highly 
impractical. In this work we develop a method of calcu- 
lation that is computationally efficient and can be used 
to study transport in large networks. Instead of finding 
the full distribution function p(f, t) we introduce a set of 
slow variables that capture the most important aspect of 
particle transport, the number of particles in each con- 
tainer A^i (t) , N2 (t) , . . . , Nm (t) , and derive equations that 
describe how they change in time. This is the central re- 
sult of the paper. 

A couple of relatedproblems have been addressed pre- 
viously in refs. 0, B 0, E3, EH Escape of a particle 
through a small hole in a cavity was studied in 0. The 
work of [8j deals with the problem of the hole connected 
to a short tube. The tube length and hole radius are 
roughly of the same size, mimicking a cell membrane 
having a thickness greater than zero. A couple of equi- 
libration cases have also been studied 0, 0, . The 
papers just indicated treat the intra container dynamics 
in much more detail than we do. In here, for simplic- 
ity reasons, the particle concentration is assumed flat in 
the container and the validity of this approximation is 
checked numerically in section rVIII In our notation, the 
studies 0, H,a,Lio, 11] can be classified as M = 1,2 cases. 
Our main interest is in networks with large M . 

This paper is organized as follows. In section ITU the 
problem is defined and the general results are stated. 
The derivation of the rate Eq. (jSJ) is explained in sec- 
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tions IIIl| and IIVI where the emptying of a single container 
into a tube, and emptying of a container into another 
container through the tube is studied. The single ex- 
ponential asymptotics of the two container system and 
related first order rate equations are found and discussed 
in section Up to this point only a two-container sys- 
tem is treated while section IVII deals with an arbitrary 
network topology. Section lYlIl elaborates on the assump- 
tion of well stirred containers. A numerical comparison to 
the diffusion equation is made. In section lVllll numerical 
case studies of various network structures are performed. 
In particular a three way junction and an example of a 
larger network are studied. The summary and outline 
of future work is given in section IIXI Technical details 
are found in the appendices. Appendix lAl describes the 
numerical procedure used for solving the rate equations. 
The rate equations for the cases studied in Section IVIIII 
are explicitly derived in Appendix [U] It can be shown 
that the presence of tube junctions can be eliminated al- 
together from the dynamical equations when the time is 
large. This is demonstrated in Appendix [D] 



II. PROBLEM DEFINITION AND MAIN 
RESULT 

Describing the particle transport in a network depicted 
in Fig. n is f ar from trivial and in order to solve the prob- 
lem several assumptions are made. We assume that (i) 
particles move solely by diffusion (the fluid in which the 
particles move stands still) and (ii) particles do not dis- 
turb each other. With these assumptions the complicated 
dynamical problem at hand is reduced to solving the time 
dependent diffusion equation: 

d tP {r,t) = V -[D(r)Vp(r,t)]. (1) 

Here p(r, t) is the concentration (particle density) and 
D(r) is the diffusion coefficient which may be position 
dependent. The walls are particle impenetrable 

d nP (r, t) = (2) 

where d n = h- V, and h is the unit vector perpendicular 
to the wall. The total number of particles is a conserved 
quantity. 

Equation JQ| could in principle be solved numerically 
using a brute force approach (e.g. the Finite Element 
Method or the Finite Difference method). However, in 
sections IIIII and IIVI we will show that it is possible to 
describe particle transport in terms of a finite number 
of variables, the number of particles in each container 
N%, ...,Nm- Also, it might be easier to understand par- 
ticle transport in such a setup. The dynamics of Ni(t) 



i = 1, . . . , M is governed by 

- e-U Cij / 4 dm (?) [a,, (t ?) 

+Kij(t-t')] 

(3) 

where 

Mi{t) = Ni{t) + N i0 S(t) (4) 

and S(t) is the Dirac delta-function, J °° dtS(t) = 1. Here 
and in the following the dot over symbol denotes time 
derivative. The connectivity matrix Cij G {0, 1} de- 
scribes how the nodes are linked (note that Ca = 0), 
dij is the radius of the tube (link) from i to j 
and Vd(Rj) is the volume of a d dimensional sphere 
Vd{r) — [27r d / 2 /dT(d/2)]r d corresponding to container j 
having radius Rj. Nio denotes Ni(t — 0). Equation © is 
derived under the assumption of ideally mixed contain- 
ers. The rate coefficients are given by 

fe=i v y / 

fe=o v y / 

£ij is the link length and is the corresponding diffu- 
sion coefficient. The theory is developed for the general 
case where the diffusion constant in each tube may be 
different. 

Figure ^ also shows the existence of tube junctions. 
They are treated by Eq. © by letting the container ra- 
dius coincide with that of the tube. Since the tubes are 
initially empty, so are junctions, N i0 = 0. Eqs. ©-© 
are the central results of this paper and their derivation 
is a major topic of the subsequent sections. 



III. EMPTYING OF A RESERVOIR THROUGH 
AN INFINITELY LONG TUBE 

To derive Eq. @ we start with the simplest possi- 
ble case and consider particle escape from a container 
through an infinitely long tube [see Fig. [21 panel (a)]. 
The main reason for this is to show how to couple the dy- 
namics of the tube and the container. Also, such setup 
captures the short time description of the full network 
problem when the particles escaping the containers do 
not yet 'feel' that the system is closed [the short time 
dynamics is contained in A(£), see Section llV] . 
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The particle concentration p(f, t) is governed by the 
diffusion equation supplemented with the boundary con- 
ditions that the walls are impenetrable and that p(f, t) 
has to vanish for x — > oo. The concentration in the 
tube and in the container are interwoven in a highly non- 
trivial way through what is occurring at the tube open- 
ing. Given that the current density out of the container 
j(0,y,z,t) is known [see Fig. [21 (a)] one could solve the 
diffusion problem and find the concentration profile in 
the container and the number of particles. Furthermore, 
one could find a relationship 



P (0-,y,z,t)=F\j{0,y,z,t)] 



(G) 



where j(0, y, z, t) = —D lim^^o- ^ ' V/?(x, y, z,t). J 7 is a 
functional that we know exists but is unlikely to be found 
in a closed analytic form. 

To find j(0; 2/> z i t) it is assumed that the concentration 
in the vicinity of the tube opening can be approximated 
with 



p(x,y,z,t) = f(y,z,t)c(x,t) x>0 



(7) 



where c(x, t) is a one-dimensional concentration and 
f(z,y,t) is a function that projects the value of c(x,t) 
onto a radial direction. Equation J7J is valid for large x 
when f(y, z, t) is constant but not in general case. Ar- 
bitrary density profile at the opening will in time smear 
out due to radial diffusion. 

By assumption, the concentration profile in the tube 
is governed by c(x,t) and it is coupled to the concen- 
tration in the container as follows. Both concentration 
and current have to be continuous as one moves from the 
container into the tube leading to 



P(0,y,z,t) = f(y,z,t)c(0,t). 



(8) 



and 



j(0,y,z,t) = f(y,z,t) lim —c(x,t) (9) 

x~*0+ OX 

Please note that lim x ^ + J^c(x,t) in Eq. 10 is a func- 
tional of c(0,t). Taking into account particle conserva- 
tion at the tube opening leads to 



dS p(0, y, z, t) — c(0, t) 



(10) 



which after using Eq. JSJ results in the condition 
J dS f(y, z,t) = 1. The problem has four unknowns 
p(0,y,z,t), j(0,y,z,t), c(0,i) and f(y,z,t) and four 
equations ©, ©-di and is fully defined. However, it 
is not tractable in this complicated form and we proceed 
to simplify it. 

Instead of p(f, t) a more useful variable is the to- 
tal number of particles in the container N(t) — 
L<o dv P(r, t) governed by 



where J(t) is the flow of particles that leave the con- 
tainer through the tube opening J(t) — J dSj(0,y,z,t). 
There are two special cases where the current J(t) can 
be determined analytically in terms of container vari- 
ables, (i) If the exit is a fully absorbing disk with ra- 
dius a, the concentration is always zero at the the in- 
terface p(0,y,z,t) — 0. In 01 it is shown that the 
current through such disk when placed at an infi- 
nite otherwise reflecting wall is Joo = AD c apoo where poo 
is the particle concentration at infinity and D c denotes 
the diffusion constant in the container. A reasonable as- 
sumption for the container, at least when the tube ra- 
dius is smaller than the radius of the container, is that 
the concentration profile far away from the exit is flat 
and can approximately be taken to be N(t)/Vd(R). Us- 
ing this for Poo yields J^t) = AD c aN{t)/V d {R). (ii) If 
the opening is completely closed, the current is zero and 
p{0,y,z,t) = N(t)/V d (R) [in such a case N(t) = N(0)]. 
A linear interpolation between (i) and (ii) yields 



p{0,y,z,t) = 



N(t) 
V d (R) 



Jit) 



Joo{t) 



(12) 



Please note that in Eq. I|12fl p(0,y,z,t) is assumed con- 
stant across the interface. This approximation is verified 
numerically in Fig. that shows f(y,z,t) ps V^_i(a) _1 . 
Also, when a <C R one has J(t)/ Joo(t) <§C 1 and second 
term in l|12f> can be neglected. This approximation is 
verified by numerical calculations in Section lYlII 

At this point, the full problem has been mapped on to a 
very simple geometry depicted in Fig. [21(b): a one dimen- 
sional line (tube) connected to a point (container). Tube 
dynamics is characterized by a one dimensional particle 
density c(x, t) along the line and all container dynamics, 
how complicated it may be, has been projected on to a 
single variable N(t). 

The only part of the problem that remains to be solved 
is the diffusion through the tube. This part of the prob- 
lem can be approximated by one dimensional diffusion 
since f(y,z,t) is constant. The constant concentration 
profile at the tube opening remains such in the tube in- 
terior (provided the tube radius does not change along 
x direction). With assumptions at hand the coupling 
Eq. JHJ) becomes 



c(0,t) = V d - 1 (a) 



N(t) 
V d (R)' 



(13) 



The concentration profile along the tube [initially empty 
c(x, 0) = 0] is given by the diffusion equation 



dc(x,t) pd 2 c(x,t) 



Of 



dx 2 



£ [0,oo) 



(14) 



supplemented with boundary conditions according to 
Eq. I|13|) and c(oo,<) = 0. The solution can be found 
by the Laplace transform method 0] and is given by 



N(t) = -J(t) 



(11) 



c(x, s) = c(0, s) e 



(15) 
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where c(x,s) — dtc(x,t)e st . Integrating Eq. © 
over the tube interface area at x = gives 



d 

Jit) = -D lim —c(x,t). 

x-^0 OX 



(16) 



Combining Eqs. JTT}, IT3|l and ltTH|l leads to a rate equa- 
tion in the Laplace transform space 



sN(s) 



N Q = - lira i VsDN(s, 

x^o V d {R) 



V d -i{a) e - Xy /I/i 



GDN(s) 



V d -i(a) 



V d (R) 

(17) 

where C[N(t)) = sN(s) — Nq. It is tempting to rewrite 
this equation in the time domain in the form of convolu- 
tion 



N(t) = - / dt'k(t')N(t-t') 



(18) 



representing a general form of a rate law, where k(t) = 
Vd-i(a) £-i jj s ^ jjowevcr, this is impossible and can 



v d (R) 

be seen in several ways. 

First, i/s has no well defined inverse Laplace transform 
and k{t) that would enter into the rate equation ((TBI) is ill- 
defined. Second, this problem could possibly be resolved 
by inverting Eq. I|15fl to obtain c(x, t) and inserting the 
result into Eq. i|16|) which leads to 



J(t) oc — lim 



d 



o dx 



dt' 



t f3/2 



e- x2/iDt ' N{t-t'). (19) 



the inverse Laplace transform 



vmg 

m 



N(t) = N exp 



Dt 



V d -i(a) 
V d {R) 



erfc 



Vd-i(o) 
V d (R) 



(20) 

Nq is the initial number of particles in the container. 



Using approximation erfc(z) 
N(t) cx ^. 



for large z gives 



Due to the complications discussed above the rate 
equation has to be stated in terms of N(t) 



N(t) 



Vrf-i(a) _ x 



V d (R) 
V d -i(a) w 



sN(s) 
dt'Af(t')A(t - 1') 



(21) 



where 



and 



V d (R) jo 
Af(t) = C- 1 [sN(s)] = N(t) + N S(t) 

A(i) = C- 1 



(22) 



(23) 



In general N(t) is unknown. To evaluate the expres- 
sion above in a way that would result in a rate equation 
involves interchanging derivation and integration. This 
is only allowed if the inte gral is uniformly convergent 
in the interval x S [0, oo) [Tj|- It is easy to see from 
Eq. I|19|) that this is not the case and the interchange is 
illegal. Another possibility is to use partial integration 
but this strategy does not work since one ends up with 
non-convergent integrals as x — > 0. Thus, Eq. I|18l) docs 
not exist for a semi-infinite case. Also, it is intuitively 
clear that one can not observe pure exponential decay 
since the system is infinite. 

For an infinite system an asymptotic rate law of the 
type N(t) oc —N(t) simply does not exist. This can also 
be seen from the exact expression for N(t) which can be 
obtained from Eq. (|17f) by sol ving for N(s) and finding where 



Please note that it is impossible to rewrite the right hand 
side of Eq. I|21|) in such a way that it would solely involve 
dependence on N(t). When the system is closed (e.g. by 
adding another container) the situation changes. 



IV. A RATE EQUATION FOR A TWO 
CONTAINER SYSTEM 

The system under consideration consists of a one di- 
mensional rod parameterized by a and £ connected to two 
ideally mixed containers having radii i?i and Ri depicted 
in Fig. H 

The diffusion equation for the tube [initially empty, 
c(x, t) = 0] connected to the two containers at x = and 
x — I is given by 

with boundary conditions analogous to Eq. (|13J) 

c (0,t) = iV 1 (t)^f, c(i,t)=N 2 (t)^$. (25) 
The solution in Laplace transform space is given by 



V d (Ri) 



V d (R 2 ) 



(26) 



sinh((x-e)y/s/D\ 



sinh x\ s 



ID) 



ssinhleWs/D) 



ssmhtty/a/D) 



Matching the current at both ends as in Eq. H16J1 

iVi(t) = D Km x ^ J^c(x,t) 
N 2 (t) = -D\im x ^ t -^c{x,t) 



(27) 



(28) 
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yields a set of rate equations in Laplace transform space 
for the two container system 



sNi(s) - Ny 



sN 2 (s) - iV 20 



-sNi(s) 



+sN 2 (s) 



-sN 2 (s 



+sNi(s) 



Vd-i(a) 




(29b) 

The terms having minus signs are the outflow while the 
ones having plus signs represent the inflow. A closer look 
at the the rate equations above reveals an important link 
to the semi-infinite case: the outflow is proportional to 
y/s for large t (small s). In other words, the semi- infinite 
case is recovered as a short time expansion of Eqs. <|29[) 
(see Section ILTTll . The physical interpretation is that ini- 
tially, the particles feel as if they were entering an in- 
finitely long tube. This implies that the outflow term 
can be divided into two parts reflecting this observation 




= A(s) + k(s). 



(30) 



A(s) is taken from Eq. I)23|l and controls the short time 
dynamics that resembles the one of the semi-infinite sys- 
tem. k(s) describes the asymptotic long time behavior 
when the particles are 'aware' of the existence of another 
side and can be found from Eq. 1|23[) and Eq. I|3U|) 



k(s) 



3 smh(lyJ7/D 



(31) 



The inflow rate is labeled a(s) 



a(s) 



smb. iyfs/D 



(32) 



The inverse Laplace transforms of A(s), k{s) and c(s) 
are shown in Eq. JSJ and depicted in Fig. [5] The long 
time behavior of n(t) and a(t) can be found from a small 
s expansion of the Eqs. (|3I|I and (|32|) 



«(*) 
a(t) 



D 

t 

D 



[l - exp (- 
Taking the limit t — > oo yields 

cr(oo) = — k(oo) 



6Dt 



)] 



(33) 



(34) 



It is more convenient to express the rate Eqs. (I29|) in 
time domain 

~ } WhJ Jo [A(t - + K(t - *' )] 



(35a) 

= f o dt'Mr{t')a{t-t') 

~wM II dt>mt,) [Ht ~ + K{t ~ t,)] 

(35b) 

n(t) and a(t) are not present in the rate equation for the 
semi-infinite case (see Eq. I21|l and arise only when the 
system is finite. 

A numerical solution to Eqs. (|35|l is shown in Fig. 
(see Appendix lAl for a more elaborate discussion regard- 
ing the numerical procedure). The number of particles 
decays exponentially which is verified in Fig.0where the 
straight line attained after some time is the evidence of 
a single exponential decay. The figure also shows a non- 
exponential regime for small t, described by A(i). Terms 
proportional to A(i) are in the following referred to as 
A-terms. 

For small t particles rush into the tube with a large 
current. At t = the current is infinite, limj^o Ni (t) — 
oo. Thus, exactly at t = it is impossible to define the 
exit rate from container and such situation extends to 
any other time instant. In strict mathematical sense, it 
is impossible to define exit rate from the container for 
any t > (when the concentration profile at the tube 
opening is different from zero). This can be illustrated 
on a simple example. Let V be a volume divided into two 
sub-volumes V% and V2 such that V\ and V2 touch each 
other and exchange particles by diffusion. The dynamical 
variables of interest are the total number of particles in 
each sub- volume Ni(t) and ^(i). The goal is to derive 
some kind of rate equation for N\ (t) and N2 (t) . We focus 
on the flow from V\ to Vi- In a small time interval e 
one will have Ni(t + e) cx iVi(t)(l - aty/e) and N 2 (t + 
e) oc N 2 (t) + /3Ni(t)\/c, where a and f3 are numerical 
constants. The effective exchange rate that describes the 
flow from V\ into V% is given by 



&2i(i) = lirn 



N 2 {t + e) - N 2 {t) 



txlime- 1/2 (36) 



which is infinite. At any time instant t an infinite amount 
of particles (per unit time) is rushing from V\ into V 2 (and 
vise versa). However, the infinite flows from V\ to V 2 and 
the other way around cancel each other out resulting in a 
finite net flow giving smooth curves for Ni(t) and N 2 (t). 
t = is special since there is no counter flow from V 2 to 
V\ which explains why lim t _>o N±(t) = 00. 

The non-exponential regime grows with tube length £ 
and might play a significant role in studying transport 
processes in networks having long connections. For large 
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times n(t) and a(t) start to dominate and one observes 
exponential decay. It will be shown in Section IVI how to 
derive rate equations that describe this regime. 



V. ANALYSIS OF THE GENERAL RATE 
EQUATION: EMERGENCE OF A SINGLE 
EXPONENTIAL SOLUTION 

It is intuitively clear that in the case of the two-node 
network discussed in previous section one should have an 
exponential decay (growth) for the number of particles in 
the container C\ (C 2 ) (see Figs. |lj]and[7J): asymptotically 
the time dependence of N\(i) and N 2 (t) is given by 



N li2 (t) = iVi, a (oo)+ A, 2 exp(-Vr) 



(37) 



where r _1 is the decay exponent that governs the late 
time asymptotics and Ai, 2 is the amplitude of decay. 
This fact is not easily predicted from the form of the 
general rate equation given in J3j. To understand the 
emergence of such behavior a more thorough investiga- 
tion of Eq. iffifll is needed. 

To obtain the exact value of the decay exponent one 
has to study the structure of poles of Ni^(s). The poles 
fully determine the form of Ni^ 2 (t) — X^Lo a p eSpt where 
a p is the residue of Ni^(s) at pole s p . The values of 
-^1.2(00) are determined by the p — term (sq = 0). 
The exponential decay rate is determined by 



T = -Si 

Rewriting Eqs. (|29[1 in matrix form yields 
sN(s) -N = MN(s) 

where 

N(s) = [N 1 (s),N 2 (s)] T , N n = [N 10 ,N 20 } T 



and 



V d (R 1 ) LOUl( i V d (R 2 )sinhq 



Vtub 



Hubc 



cothq 



(38) 

(39) 
(40) 

(41) 



with q 2 = s£ 2 /D and V tuhc = V d -i(a)£. The 
poles are calculated from det(s — M) = since 
N(s) = (s - M^Nq and (s - M)' 1 cx 1/ det(s - M). 
Evaluating det(s — M) = gives 



q 2 + q 



Vtu 



be 



be 



V d (R!) V d {R 2 ) 



coth q 



V 2 

tube q 



(42) 

Equation 142|l is a transcendental equation and has 
many solutions q p that determine the value of the poles 
s p = q 2 D / i 2 where p — 1, 2, . . . , oo and in particular 



si 



q\D 

I 2 



(43) 



that, together with Eqs. (|38|) . gives a relationship be- 
tween qi and the decay rate r _1 = —q\Djl 2 . Note 
that q 2 has been factored out from Eq. (|42|l and sq = 
(qo = 0) is the additional pole that determines the values 
of -ZVi^oo). Also, note that q\ depends only on parame- 
ters describing geometry of the network. 

Apart from determining the structure of the poles, 
Eqs. I|39l) - (|41|l are a good starting point for classifying 
various schemes for obtaining approximative forms of 
the rate Eqs. (|35|l . Equations (|35|l do not have the 
form of the general rate law stated in Eq. (|18fl (due to 
the presence of the A-terms). Such a rate law would 
be easier to understand intuitively. For example, the 
emergence of the single exponential decay could be 
seen more easily in Eq. (|18fl than in Eqs. I|35|l . Also, 
an approximative form might be easier to implement 
numerically, though at the cost of a lower accuracy at 
the end. The idea is to perform a small s expansion 
of Eq. based on a desired accuracy. In here we 

consider two cases. 

(a) Lowest order expansion: Performing the expan- 



q coth q w 1 



sinh q 



(44) 



of M. in Eq. 14111 leads to a matrix that is constant, and 
taking the inverse Laplace transform of Eq. 1)39(1 gives the 
following set of rate equations 



N 1 (t) = -N 2 (t) =VtubeS 



N 2 (t) N^t) 



V d {R 2 ) V d (Ri) 



(45) 



These equations were already stated in ref. It is in- 
teresting to see that they emerge as a special case of the 
scheme discussed here. Also, Eq. I(45|l can be obtained by 
following another route. Performing partial integration 
of Eqs. (|35() with terms containing A(t) omitted leads to 
the exactly the same form of rate equations as given in 
(|45H . This procedure is discussed below. 

The A-terms are only present when the system is in- 
finite. Since the problem is finite, terms proportional to 
A(t) will be sub-leading for large t. This can be seen 
from a small s (large t) expansion of Eqs. (|23|) and l|31|) 
j3t|. Also, partial integration of the A-terms is impossi- 
ble since the derivative of A(t) is proportional to i~ 3 / 2 
and diverges when t — > 0. 

Omitting A-terms in Eqs. I|35|l. and performing partial 
integration leads to 



Ni(t) = 



JVa(t) 



(46a) 



(46b) 
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Since &(t) is peaked for small t (see Fig. 01, the contri- 
bution to the integrals (convolutions) stems mainly from 
small values of t' which justifies the following approxima- 
tion 



dt'N 1<2 {t-t')&{t') m JVi, 2 (i) / dt'a(t')=N h2 (t)a(t). 
o Jo 

(47) 

where <r(0) = was used. The same applies for n{t). 
Using Eq. gTJ in g?| leads to 

Mt) = ^#T^(*)^(*)-^#T^(*)«(*) (48a) 



V d (R 2 



= Hy^^^li^^^ (48b) 

Inserting k(oo) and cr(oo) found in Eq. Il.'i 111 into i|48|) 
yields l|45|l . This example shows how the A-terms dis- 
appears from the description when the system is finite. 
However, contrary to the partial integration method, the 
expansion of Eq. (|39[1 gives a more systematic and con- 
trolled approach. 

Equation (|45|) is simple and computationally efficient. 
It could be easily used to describe large networks. How- 
ever, it has several drawbacks that can be identified. The 
solution to Eq. 145|) is given by 

N h2 (t) ~ JVi, a (oo) = [JVi, a (0) - ^1,2(00)] exp (-t/n,«) ■ 

(49) 

The decay rate r a 1 = —q 2 a D/£ 2 is determined by 



Ql,a 



V tuhe [V d (Rl)+V d (R 2 )} 



V d (Ri)V d (R 2 ) 
The number of particles in each container as t — > 00 is 



(50) 



N 



10 



20 



JVi(oo) _ N 2 (oo) _ 
V d (Ri) ~ V d (R 2 ) ~ V d (Ri) + V d (R 2 )' 



(51) 



Equation Q51JI is not correct. The correct values for 
-^1,2(00) are given by 



N 



10 



211 



iVi(oo) _ jV 2 (oo) _ 

Vd(Ri) ~ V d {R 2 ) ~ Vd(Ri) + V d (R 2 ) + U tubc 



(52) 



The discrepancies between Eq. (fSTJ) and (p)2*|) become 
increasingly important for long tubes which are likely to 
occur in large networks. For example, in a case where 
the tube and reservoir volumes are equal, Eq. I|51|) 
predicts A/i(oo) = N 2 (oo) = N tot /2, N tot = N 10 + N 20 , 
while the exact result from Eq. I|52|l is Ntot/3- The 
particle decay exponent given in Eq. (|50|l only holds 
when T4ubc — > 0. It strongly deviates from the exact 
value when the tube is long (see Fig. [SJ . 

(b) Higher order expansion: Using the expansion 

(? coth« ? «l + £ -^«1 (53) 



for inserting in Eq. (|39fl and taking inverse Laplace 
transform, leads to a set of rate equations (given in Ap- 
pendix [BJ that are unsatisfactory due to the following 



reasons. First, they predict a spurious jump in Ni^ 2 (t) 
as t — > 0, and the limiting values N\^ 2 {oo) are not cor- 
rect. Second, the rate exponent that results from these 
equations is not that accurate. This particular example 
shows that it is important to have a balanced expansion 
for elements of M.. For example, instead of expanding 
q coth q directly one has to expand sinh q and cosh q sep- 
arately in such a way that same powers in the nominator 
and denominator are obtained. When this strategy is fol- 
lowed a much better approximation is obtained as shown 
bellow. 

The next order expansion, gives correct limits for 
Ni^(t) when t — > and t — ► 00 and leads to a rela- 
tively accurate value for the decay exponent. Using the 
expansion 



n cnth a ~ 1+q2 J 2 1 ~ * 

gcotng- 1+q i/ 6 S!^ ~ i+q 2 /e 
in M. gives the following set of equations 
3D Kubc 



(54) 



Ni(t) = -JVi(t)- 



P V d (Ri 



+- 



12D 2 V tub 



Na{t) 



£ 4 V d (R 

6D 2 Kube 

" ** V d (R 2 

3D Vtub. 



1) Jo 

f dt'N 2 (t-t')exp 
Jo 



6Dt' 
6Dt' 



(55a) 



12D 2 Vtu 



e 4 v d {R 2 ) 

6D 2 U tub 



P V d {R 2 ) 

'" / dt'N 2 {t-t')exp 
Jo 



V d (R 



\ [ dt'Ni(t-t')exp 
1) Jo 



6Dt' 
QDt' 



(55b) 



Solving det(s — A4) = to get hold of the decay ex- 
ponent in analytical form becomes in this case rather te- 
dious since finding the value q\ amounts to finding a root 
of a fourth degree polynomial. The calculation simpli- 
fies somewhat if equal container volumes are considered, 
Vd(Ri) = Vd(R 2 ) = Vd- In such a case one has 



lib = -W7[3(2U d + U tubc 



v/3(12^ - W d V tuhc + 3V 2 uhc ) 



(56) 



The main findings of this section are summarized in 
Figs. H3E|and|Hl Figures H3 and depict a numerical so- 
lution of Eq. (|35|l (solid line) compared with the approx- 
imations discussed in this section for a case where the 
tube and container volumes are equal. Figure |H1 shows a 
detailed analysis of the decay rate. 

For very short times there is a difference between 
Eqs. H55fl and (|35|l in Fig. and [7| These arise due to 
the partial elimination of A-terms. For example, Eq. I|55|) 
does not predict lim t _>o Ni (t) = 00 [Fig. El the dashed 
line lies above the solid line near t — for curves depict- 
ing Nx(t)]. This is the reason why the curve for Ni(t) 
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obtained from Eq. I|55|l underestimates the emptying of 
container C%. This effect is more pronounced for the 
Ni(t) coming from Eq. (|49(1 . There, the A-terms are 
eliminated altogether [Fig. HJ1 the dotted line depicting 
Ni(t) lies above solid and dashed lines]. Also, Figure 
shows that there is a large error in ./Vi^oo) for curves 
obtained by Eqs. (gSJ- 

In Fig.dthe natural logarithm of [Ni(t)-Ni(oo)]/N tot 
is shown. For short times the dynamics is not exponen- 
tial but after some time a straight line is attained which 
is the evidence of single exponential behavior. The curve 
corresponding to Eq. I|49l) (dotted line) does not predict 
the correct decay exponent which is manifested in a dif- 
ferent slope. The decay exponent predicted by Eq. 
is a better estimate for the decay rate: the slopes of the 
dashed and solid lines more or less coincide. This fact is 
shown more clearly in Fig. [H] 

Figure 03 depicts the dependence of q\ as a function 
of the tube volume. The numerical solution to Eq. i|42|) 
(solid line), which gives the exact value of the decay ex- 
ponent, is compared to the values of q\ a (dotted line) 
and q\ b (dashed line). All three cases work well when 
Vt u be — >■ 0. As the tube volume increases q\ a deviates 
more and more from q\ . The same holds for q\ b , though 
its value lies much closer to q\. For example, when all 
volumes are equal q\ = —1.71 and q\ h = —1.63 while 

Via = -2. 

In this section methods of finding rate equations and 
decay exponents for the two container problem was de- 
duced. Figures |H1 13 and show that for increasing tube 
volumes the rate equations given in Eq. I|49l) are not ca- 
pable of describing the dynamics. The approximation 
given in Eq. (|55J) works better. For short time dynam- 
ics none of the developed methods are valid and the full 
rate Eq. @ is the only alternative. In the subsequent 
section all rate equations discussed up to this point will 
be extended to work for any network structure. 



VI. THE GENERAL EXPRESSION 



This implies 

N t (t) = &j [ IN i-< - OUTj^] , i = 1, . . . , M 

(58) 

where is the conductivity matrix discussed in Sec- 
tion [n] The final result is stated in Eq. ©. 

It was argued earlier that if the tube volume is small 
a very simple first order rate equation can be stated. It 
is found in Eq. (|45|l . Extending this equation to the case 
of arbitrary topology yields 



Nj(t) Ni(t) 



V d (Rj) V d (Ri 



Note the symmetry relations 



(59) 



and 



Dij — Dji. 



If a more sensitive solution is desired, a set of equa- 
tions of the type ifSTj) is suggested. An extension of this 
equation is shown below 



OUT^ = Ni(t) 



3Djj Vd-ijctij) 
tii V d (Ri) 



' dt exp ( — ^ — )Ni(t-t) 



1% V d (R t ) , 



I 2 - 



(60a) 



(60b) 



Inserting Eq. (|60|l in Eq. (|58|l yields an approximative 
form of rate equations for an arbitrary network. 

In previous section we investigated differences between 
Eqs. (|57|) . (|59l and Eq. JBDJl using the two- node network 
as a study case (M = 2). It is expected that the findings 
of previous section also apply for larger networks. This 
analysis is not conducted here. Having general expres- 
sions at hand more complicated network structures can 
be investigated. However, the assumption of well stirred 
containers will be discussed first since it is crucial for the 
derivation of the rate equations P|l. 



In this section results and methods obtained and de- 
veloped for two-node network will be extended to work 
for any network structure. The complete dynamics for 
the two-node network is formulated in Eq. i|35|) . For an 
arbitrary network the outflow (OUT) from container i to 
container j is proportional to Ay (t) + Kij (t) and the in- 
flow (IN) from container j to container i is proportional 
to aji(t): 



01 T , = Vd ~ l{ ^ ) I dt'Mi{t')[^{t-t') 



IN,- 



V d (Ri) j 
+K ij (t-t')] (57a) 

%M f dt'M 3 it')°At-t'). (57b) 

Vd{Kj) Jo 



VII. THE ASSUMPTION OF IDEALLY MIXED 
CONTAINERS 

An ideally mixed container has no concentration gradi- 
ents. If a particle has entered, it can be found anywhere 
within the compartment with equal probability. In real- 
ity, a diffusing particle examines the compartment in a 
random walk fashion until an opening is found. There 
it has a possibility to escape and change the concentra- 
tion. If the tube radius is small (o/fl«l),a significantly 
longer time is required to find an opening and escape than 
to examine the majority of the compartment. The time 
needed to examine the majority of the compartment is 
called mixing time and is given by r m ; x = Sfj. The time 
of finding a specific place or target having radius a is 
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given by r tar g t = • For cubic compartments the 

radius of the sphere R is replaced by the edge length, in 
this case 2R (see Fig.^J). Thus, we argue that if a <C R 
then T m i x <C Te SC ape and the containers can be consid- 
ered ideally mixed at all times (^3 ■ This is supported by 
numerics. 

Figure 03 (b) shows a numerical solution of the diffu- 
sion equation in two dimensions in a geometry depicted 
in Fig. [5] (a). The solution clearly shows that the as- 
sumption of well stirred containers becomes very good 
for a/R <C 1, even for skewed initial distributions. The 
solution to the diffusion equation was found using a stan- 
dard implicit finite difference discretization method |16| . 

In Fig. El (b) one observes a systematic discrepancy: 
the assumption of ideally mixed containers tends to over- 
estimate the decay rate. This derives from the fact that 
the assumption of well stirred containers over estimates 
the number of particles at the tube inlet leading to a 
larger exit rate. 



VIII. CASE STUDIES 

Up to this point the two node network was the 
only example discussed. It was used as elementary 
building block when constructing the rate equation Jjy|. 
Such network is rather simple and does not offer any 
spectacular behavior. In this section more complicated 
structures will be studied that are shown in Figs. El 
and 1171 Two realizations of a three-node network are 
studied first, case 1 and 2, depicted in Fig. El (a) and 
(b). They possess one more level of complexity than the 
two- node network shown in Fig. 0] Case 3, Fig. El ( c )> 
is a four-node network that has a T-shape structure. 
Case 4, Fig. E3(d), has the shape of a star and involves 
a tube junction. It will be shown that the transport 
properties of this structure can be controlled in such 
way that it might serve as a diffusion-based transistor. 
These three and four node- networks, despite being 
rather simple, exhibit a large variety of outputs. Case 5, 
Fig. El demonstrates the transport properties for larger 
networks that are impossible to predict without using 
Eq. ©. When solving Eq. numerically the con- 
tainers and the tubes are considered three dimensional. 
Moreover, the diffusion coefficients Dij and the tube 
radii were kept the same for all tubes. All graphs are 
scaled with the total number of particles in the system 

Case 1: The Line. Three reservoirs are lined up on 
a straight line as depicted in Fig. El ( a )- The transport 
equations for this system can easily be written down us- 
ing Eq. © and are listed in Appendix IC II A numerical 
solution is shown in Fig. ^2 This three-node system ex- 
hibits a new characteristic that can not be found in the 
two node network. The curve depicting the number of 
particles in the middle container (solid line) has a maxi- 
mum (an extremum) point. 



The extremum point is a manifestation of an un- 
balance between inflow and outflow in the middle 
container. This unbalance derives from an asymmetry 
in the structure and arise only when £12 < £23 or 
Vi < V3. This provides a possibility of designing the 
output pattern through simple geometrical changes in 
the structure. Fig. 1121 is a simple demonstration of this 
design possibility where £12 is varied so that the arrival 
time of the maximum of N^it), depicting the number 
of particles in the middle compartment C2, is changed. 
The picture also shows that an increase in £12 results in 
both an increased arrival time and a wider peak. The 
height of the maximum is controlled by the value V2. A 
decrease in V2 suppresses the peak and vice versa. 

Case 2: The Triangle. The connectivity of the line 
is changed by adding an extra link between containers 
C\ and C3 so that it forms the shape of a triangle, as 
shown in Fig. I1UI (b). The rate equations are derived 
in Appendix IC 21 and a numerical solution is shown in 
Fig. El The initial distribution of particles, see the inset 
in Fig. 1131 is chosen in such a way that a minimum is 
produced in the curve depicting the number of particles 
in container C2. 

The extremum point can be enhanced or removed com- 
pletely in the same way as was demonstrated for case 1. 
Such minimum will be absent unless the geometry and 
initial distribution of particles are tailored in a specific 
way. In general, such sensitivity of geometrical changes 
and changes in the location where particles are injected 
is observed for all cases studied in this section. 

The triangle structure studied here exhibits shorter 
equilibrium process than the linear structure considered 
previously [see Fig. El (a)]- In the case of a triangle 
structure particles can spread more efficiently due to the 
additional routing possibility C\ — > C3. 

Case 3: The T-Network. It is interesting to see how 
the behavior of the cases 1 and 2 changes when a new 
node is added to the network. In here we study a situ- 
ation where an extra node C4 is connected to container 
C2 in the structure depicted in Fig. El ( a )- In sucn a way 
one gets a T-shaped (star) network shown in Fig. 1101(c). 
This alternation of structure leads to a significant change 
in behavior, as shown in Fig. El When compared to the 
cases 1 (one maximum) and 2 (one minimum) the curve 
depicting iV 2 (i) exhibits an additional extremum point: 
both minimum and maximum are present simultaneously. 
The right inset in Fig. 1141 emphasizes this fact. 

This scheme could be carried out further adding on 
more and more reservoirs and adjusting the lengths and 
the initial distribution so that the peaks arrive in con- 
secutive order, possibly produce an wave-like behavior. 
However, since the spread of the peaks increases with 
increasing tube length, it might be numerically quite 
difficult to see when the the extremum points occur or 
even if they actually exist. 
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Case 4: The Junction. The next interesting network 
to consider is the one with a junction present as depicted 
in Fig. EH (d). A particular example of a three way junc- 
tion is studied. The network is built up by three reser- 
voirs and three tubes. The ends of the tubes coincide 
to form a junction. To obtain the transport properties 
of such a network we start from the structure studied in 
case 3 shown in Fig. ^J|(c). The junction is obtained by 
reducing the radius of container Ci in the middle until 
its radius is roughly equal to the radii of the surround- 
ing tubes, see Fig. ED The rate equations describing 
the junction properties have the same form as the equa- 
tions describing case 3 (listed in the appendix IC 31) with 
the substitution Vd{R-i) ~ > Vd(a). The equations are not 
given in order to save space. 

A numerical solution of the rate equations for the junc- 
tion is shown in Fig. ^| This structure allows control 
of the particle flow between V\ and V3 by adjusting the 
volume V2 and the length £24 (see Fig. 11511 . This setup 
could function as a diffusion based transistor. For exam- 
ple, by making £24 shorter than £34, compartment C2 will 
initially attract diffusing particles from C\ to a greater 
extent than C3 causing a time delay in the particle arrival 
into container C3. 

This fact is illustrated in Fig. 1161 where the maximum 
in the curve for the number of particles in container C2 
(solid line) indicates an initial accumulation of particles 
in C2' N2(t)/Ntot rises from to 0.5 in the interval t = 
to Dt/£ 2 — 5. In this interval container C2 accumulates 
the majority of the particles released from container C\. 
After the peak has been reached the particles accumu- 
lated in C2 are released into C3 : N2 it) /N to t continuously 
drops from value of 0.5 after Dt/£ 2 = 5. By changing £24 
the curve for A^(t) can be manipulated exactly in the 
same way as done in Fig. 1121 but such analysis is not 
repeated. 

For large networks involving junctions it might be de- 
sirable to decrease the number of equations required to 
solve the transport problem. It is demonstrated in Ap- 
pendix[D]that the presence of the junctions can be elim- 
inated altogether when investigating dynamics in the 
t — *■ 00 regime for structures having large container vol- 
umes (a < R4, i = 1, ...,M). Equations (|D8)l - IID10|) 
show this explicitly for the three node junction studied 
here. 

The four cases studied up to now show that the trans- 
port dynamics is very sensitive to geometrical changes 
and to the locations where particles are injected. Small 
variations in tube lengths and container volumes lead 
to unpredictable changes in curves depicting the time 
dependence of the number of particles in each container 
(Figs. llllll3lini and [T^|) . Also, the shapes of the curves 
differs significantly from a single exponential decay . 

Case 5: Large Network. Following the methods de- 
scribed in this paper one could easily study diffusive 
transport in networks containing hundreds or thousands 
of containers, tubes and junctions. The computational 



cost scales linearly with both the number of containers 
and number of tubes (assuming full connectivity of the 
containers). We do not show explicit example with such 
large number of containers. To demonstrate the power 
of the method we study a more pedagogical example, the 
case of a network that is built up by seven containers, 
9 tubes and a 4-way junction. In contrast to the previ- 
ous cases shown in Fig. ^| it is impossible to predict the 
transport behavior of such a network without a numerical 
calculation. 

Figure ITS1 panels (a)-(c), shows numerical solutions 
of Eq. © for the network depicted in Fig. El Different 
initial distributions are used and are introduced into inset 
of all figures. The darker the container appears the more 
particles are injected into it. The dynamics is evidently 
quite complex and all possible characteristics that were 
forced upon the other cases are present. There is an 
exponential growth and decay as well as curves having 
one or more extremum points. The different transport 
behavior shown in Fig. ll8l stem only from different initial 
conditions. If the structure no longer remained fixed even 
more complicated patterns could be produced, only the 
imagination sets the limits. 



IX. CONCLUDING REMARKS 

We introduced a generic model for a diffusive particle 
transport in large networks made of containers and tubes. 
The diffusion equation that describes the distribution of 
particles p(f, t) throughout the network was taken as a 
starting point. Instead of calculating p(r, t) explicitly, we 
followed another route and developed a theoretical tech- 
nique to solve the transport problem using finite number 
of variables Ni(i), . . . , NM(t) that describe the number 
of particles in each container. First, a set of rate equa- 
tions was derived for the two-node network and second, 
they were generalized to work for an arbitrary network 
structure. In such a way we obtained the rate equations 
that govern dynamics of Nx (t), . . . , Nm (t) ■ These equa- 
tions are summarized in Eq. J2J) and are the central result 
of the paper. 

The transport equations were found by study the ex- 
change of chemicals between the container and the tube. 
It was demonstrated in Section ITTll how to couple the dy- 
namics in the containers and tubes, as stated in Eqs. © 
and (|51)- (|10f) . However, the coupling is too complicated 
to be carried out in practice in the original form. Several 
approximations were made in order to make such scheme 
doable. 

The tubes were assumed to be one dimensional lines 
(see Eq. |SJ, and the transport in the tubes was described 
in terms of a one dimensional diffusion problem involv- 
ing the concentration profile along the line c(x,t). The 
expression for c(x, t) can be found analytically using e.g. 
Laplace transform technique. In such a way the tubes 
were eliminated from the problem. 

We have considered diffusive non-interacting point par- 
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tides which is a plausible assumption for dilute solutions. 
One could consider the situation when the particles dis- 
turb each other. In that respect, the region of the tube 
interior is the most critical since the tubes can be very 
narrow and exclusion effects will be mostly pronounced 
in there. Such effects could be added into the theo- 
retical description by using results obtained from stud- 
ies on diffusion with exclusion in one dimensional sys- 
tems 0, 0, . We expect a very different transport 
behavior when the diameters of the tubes becomes com- 
parable in size to that of the particles. 

The dynamics in the container is not tractable ana- 
lytically but it was argued that when the tube radii are 
smaller than any other length scales in the system (e.g. 
tube lengths or container radii) the containers can be 
treated as ideally mixed at all times. This assumption 
was verified numerically in Section fVIII and simplifies all 
intra-container dynamics to one dynamical variable: the 
total number of particles in a given compartment. 

Evidently much of the container dynamics has been 
neglected but the coupling is formulated in such a way 
that a more detailed description can be developed should 
there be a need for that. For example, the container dy- 
namics could be better treated by using the techniques 
presented in e.g. 9], or by further exploration of the cou- 
pling equations |J?J| and l|S|l- i|l(J|) . For example, one could 
keep the second term in the right hand side of Eq. i|12|) 
and use p(0,y,z,t) = N{t)/V d {R) - J(t)/4D c a instead 
of p(0,y, z,t) = N(t)/Vd(R). This procedure would lead 
to similar rate equations as presented here with different 
forms for A(i), n(t) and a(t). 

Initially the particles feel as if they are escaping from 
the container into an infinitely long tube. In this regime 
the number of particles in the container decays non- 
exponentially. We have identified terms in the rate equa- 
tions describing this behavior: all terms in Eq. JJJl pro- 
portional to A(t) dominate when time is small. This non- 
exponential regime grows with increasing tube length and 
crossover time for this regime scales as £ 2 /D. Also, terms 
containing A(t) contain solely dependence on Ni(t) and 
can not be rewritten in the form that would involve Ni(t). 
Accordingly, it it impossible to rewrite Eq. (j3Jl so that it 
adopts a form of a general rate law (see Eq. I18fl . These 
issues were discussed in Section ITTT1 The bottle neck lies 
in the definition of transport rate which, in principle, is 
an ill-defined quantity (see discussions at the end of sec- 
tions HUl and H2|) ■ 

Eventually, for large times the decay is exponential. In 
such regime term proportional to A(t) can be neglected 
in the rate Eq. @. Other terms proportional to n(t) and 
a(t) can be rewritten in the form of a general rate law, 
leading to Eqs. (O and (JHOJ). We showed that Eq. is 
a special case of the general approximation scheme devel- 
oped in Section starting from rate Eq. J3J in Laplace 
transform space we developed a series of approximations 
that can be used to systematically describe the asymp- 
totic regime, resulting in Eq. IjtiUfl . Also, we have devel- 
oped a procedure that can be used to eliminate junction 



points for large times which reduces the number of vari- 
ables further (see Appendix |D"|) . 

Already simple case studies that were used to illustrate 
the workings of the method exhibit interesting behavior. 
For example, one can identify three types of curves that 
appear in the plots depicting the time dependence of the 
particle number in each container (Figs. 1111-1141 and llGI) . 
Type I curves occur for the two node network. These lack 
extremum points and the particle number either strictly 
rises or drops to saturate to asymptotic values. Type 
II curves have one maximum or minimum, and type III 
curves can have more and these are the most interesting. 
Type II and type III curves normally describe the particle 
number for the container in the network interior. 

The existence of type I curves for large networks sug- 
gests that it could be possible to understand the trans- 
port between two nodes in terms of an effective two node 
network where a complicated structure of links and con- 
tainers in between two nodes is mapped onto an effec- 
tive link connecting them. One can raise a more general 
question: what is the smallest network that would have 
the same transport properties as some sub-structure of 
a given large network? For example, if one is interested 
in only three nodes of the structure depicted in Fig. El 
e.g. C5, C7 and Cg, is there a star-type network, e.g. as 
the one depicted in Fig. ^| (c) or (d), that would have 
equivalent transport properties? 

The presence of type II and type III curves indicates 
the possibility that there might be curves that posses a 
larger number of extremum points. These are likely to 
occur in larger networks. We can of course manually en- 
hance certain properties as height and width of peaks. 
This will however become more and more complicated 
as the network size increases. Problem is that the peaks 
that occur later have larger width and might be harder 
to see. For design purposes, it is therefore necessary to 
build a learning mechanism or a search engine, on top of 
our existing software, to select certain characteristics in 
the curves depicting time evolution of N%, Nm- Fur- 
thermore, to exploit such effects one has to amplify them 
in some way. At the moment our study deals with trans- 
port only, but reactions in the containers can be included 
as well, and they could be tailored to amplify such effects 
(e.g. by choosing a reaction of enzymatic type). 

The techniques developed in this work can be used to 
study transport in various systems. In the following we 
give a few examples. 

(i) Although we exclusively study transport, our model 
can serve as a platform for reaction-diffusion-based bio- 
computing devices Jl 11 IIEI El El El El 
El El El El El EJEH1 In particular, our work 
is applicable to studies of the reaction-diffusion neuron 
El Hi El El El E! El- The reaction-diffusion neuron 
is an 2D array of compartments that exchange chemicals 
by diffusion. 

(ii) A large number of processes happening in the cell 
are governed by transport of reactants and chemical re- 
actions. In order to avoid a need for excessive storage 
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facilities the chemical compounds are routed in an or- 
derly fashion between various places within and between 
the cells and the chemical components arrive exactly at 
the right place at the right time [TH l3f| . The setup in 
Fig. captures this aspect of the cell interior. 

(iii) The transport on abstract mathematical networks 
(nodes and links) has been studied extensively [36l \3l\ . 
Such studies are geometry free with emphasis on the 
topology of the network graph (the connectivity patter, 
the average number of neighbors etc.). The techniques 
developed in this work could be used to account for the 
fact that the links between the nodes have physical length 
and the transport along the links is not instantaneous. 

In summary, the work we have presented is a step to- 
wards understanding the transport properties of large 
networks where geometrical concepts such as the length 
of the tubes play an important role. The setup employed 
in this study is rather simple. In order to be able to 
focus on issues related to transport, the reactions are 
totally omitted. The concentration profile in contain- 
ers is assumed flat and this is good approximation when 
tubes are thin. Already simple examples of networks we 
studied show a number of interesting properties. For ex- 
ample, transport properties of the networks exhibit large 
sensitivity to the geometrical changes in the structure. 
Also, one can adjust structure to obtain wave-like be- 
havior (with one or two extremum points) in the curves 
that depict number of particles in containers. When 
the complexity of the network increases one can expect 
even more complicated behavior with larger number of 
extrema points. The setup we use is generic and it is 
possible to expand the model in many ways, e.g. by im- 
proving description of intra-container dynamics, incor- 
porating reactions, or allowing particles to disturb each 
other. It will be an interesting problem to try to explore 
these questions further. 
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APPENDIX A: NUMERICAL CONSIDERATIONS 

In this section the numerical solution of Eq. © is dis- 
cussed. In Section ITTT1 and Section llVl it was shown that 
the time derivatives Ni(t), . . . ,NM(t) can be infinite at 
t = which may cause numerical difficulties. However, 
the singular part of the derivative can be factored out by 
making the substitution 



i = 1, 



(Al) 



where ipi (t) is a smooth function of time. For small t the 
particles in the containers do not (yet) 'feel' the presence 



of another side and behave as if entering an infinitely 
long tube. An expression describing the transport be- 
havior for for such a case was found analytically and is 
stated in Eq. JSUJ. The time derivative of Eq. 1(2*0)1 is 
proportional to i -1 / 2 for small t. Inserting Eq. (|A1I) in 
Eq. © yields an equation for ipi(t) that was used for 
numerical calculations 



c 



Vd-i(aij) 
V d (R t ) 



|/o|^(t')«ii(«-*') 



(A2) 

Note that the expression for Ay(i) has been inserted into 
the integrals. Combining Eq. i|20[l and (|A1|> leads to 

Mt) = N l0 - (Trt^Niit), i = 1, . . . , M. (A3) 

and sets the initial condition ipi(0) = Nm. 

From Eq. (|A2|1 two types of integrals can be identified: 



4 n bP] - 



(A4) 



where tp(t) is non-singular in the range of integration. 
The quadrature formulas derived to solve l\ n \tp] and 
J2" [ip] are based on the methods described in [38j |. The 
idea is that the singular part of the integrand, 2 7-1 / 2 and 
[t'(t n — t')]^ 1 / 2 respectively, is treated exactly while the 
smooth part ipi(t') is linearly interpolated between ti and 
tj+i. The resulting quadrature formulas are of the form 



it M 



1,2 



(A5) 



where 



are weights, n — 1,2... and t n = nh. This 



quadrature formula becomes exact when ij)(t) is piecewise 
linear. Calculating weights for Ii[4>] yields 



a0 = yjn - 1 - (n - 2) 



arcsin 1 / — 



- ~i + 1) + y/U + -i - 1) 



-Vi(n-J') + 2 (i + 1 -f) 



—2 arcsin \/4 + arcsin \ ^ii 

y n y n 

s/n — 1 + 7r (l — + (n — 2) arctan \/n — 1. 



(A6) 
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Calculating weights for I 2 [ip] yields 



WnO 



U'r, 



3 

Ah 1 ' 2 



[(j + I) 3 / 2 - 2// 2 + (j - I) 3 / 2 ] 



Wnn = 2hVl [„3/a _ 3(n _ !)i/2 + 2(n _ 1)3 /2] _ 

Finally, an expression for the total number of particles 
is found by integrating Eq. I|A1[) 

r f dt' 

Ni(t) =Ni(Q)- / — Vi(t'), i = l,...,Af. (A8) 
Jo 

Ni(t), . . . , Nj\i(t) are found by using the quadrature for- 
mula derived for [-0] . 



APPENDIX C: RATE EQUATIONS FOR THE 
CASE STUDIES 



This appendix explains in a more detail how to de- 
rive the rate equations used in the cases studies in Sec- 
tion lVIlfl The equations arc obtained from Eq. © . Also, 
Eq. (|59"j) is used to illustrate the impacts on dynamics 
from changes in the network structure in a less compli- 
cated form. The Afi(t), Aij(t), K%j{t) and <Jij{t) are de- 
fined in Section [H] The initial distribution is set to be 
Nj(0) = Nj where j = 1,...,M and M is the total 
number of nodes in the system. Note the symmetry re- 
lations ay = dji, £{j = £ji and Dij = Dji which implies 
that «ij(t) = Kji(t), Aij(t) = Aji(t) and a i3 (t) = (Jji{t) 
(see Eq. EJ. 



APPENDIX B: FINDING A SET OF 
APPROXIMATIVE RATE EQUATIONS FROM 
AN EXPANSION OF EQ. IN THE 
VARIABLE s. 

In this section, details for obtaining a set of a first order 
rate equations from Eq. I|39|l will be outlined. The dy- 
namical equations are found by approximating M. given 
in Eq. (|41|l . M is approximated by using the expan- 
sion stated in Eq. I|53|) . The inverse Laplace transform of 
Eq. H39|) after approximation reads 



N a (t) = 



SVd(Ri)Vtu bc D 

3Vd(fli) + y tube p 
-S{t) Kubc 



N 2 (t) JVi(t) 



3V d (Ri) + V tvS 



V d {R 2 ) 
-N w 



V d {Ri 



3V d (R 2 )Vtubc D 

3Vd(Ra) + Kube P 

-5(t)- Vtuhc 



m(t) N 2 (t) 



V d (Ri) 
N 20 . 



V d (R 2 



(Bla) 



(Bib) 



2 

01,5' 



3V d (R 2 ) + Kubc 

The decay rate predicted by these equations is given by 
rfj, = —£ 2 /qf b ,D where 

VtuUVd(Ri) + Vd(R*)] (B) 

V d {R 2 )[V d {R x ) + V tuhe /2>] 1 ' 

This decay exponent is not adequate. For example, in 
the case where all volumes are equal q\ h , = —1.5 while 
the exact value is <? xact = —1-71. 

The rate Eq. (|B1|I can not describe the behavior of 
N\ t 2{t) as t — ► 00 and t — > 0: the values for ^1,2(00) are 
given by 

ATi(oo) _ N 2 (oo) _ N 10 + N 20 



V d (Ri) VdiRi) V d {R x ) + V d {R 2 )- 
and for iV 12 (0 + ) there is a sudden jump 



\Vtu 



be 



^(0 H 



1 



1,2 



(B3) 



(B4) 



V tuhe /W d (Ri) 

and Ni(0 + ) ^ Ni(0). This is not satisfactory and another 
type of expansion is needed if correct limits and better 
decay rates are to be found. 



The Line 



The structure of the linear network studied here is 
shown in Fig. ^] (a). Containers C\ and C3 are con- 
nected only to the middle container C 2 . Therefore, the 
rate equations describing the dynamics in each container, 
N±(t) and N 3 (t) respectively, have the similar form: 



Ni(t) 



+A i2 {t-t')], i = 1,3. 



(CI) 



The middle container C 2 is connected both to container 
C\ and C3. The rate equations for N 2 (t) reads 



N 2 (t) 



Y^^ft dm{ t>)a 12 (t-t') 



+ 



Vd-l("32) ft A+ l 

V d (R 3 ) Jo 



fodt'Af 3 (t')a 32 (t-t') 



- Y WM 1 ft *') + A 2i(* - t')] 

- Y Tim 1 Jo dt'M 2 {t')[K 23 {t t') + A 23 (t t')}. 

(C2) 

Equations (|C1(I and (|C2|I have a quite complicated struc- 
ture and impacts on the dynamics of N\(t), N 2 (t) and 
N 3 (t) due to changes in the structure (e.g tube lengths 
and container volumes) are not easily predicted. Simpler 
expressions can be obtained by using Eq. (|59|) which is 
valid in the large time limit under the assumption that 
the number of particles in the tubes is negligible over 
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time. Applying Eq. leads to 



Ni(t) 



N 2 (t) 



V d (R 2 ) In 



N 2 (t) 



Vd-l("i2) D N (,\ 



1,3 



Vd-i(ai2) D 

V d (R 2 ) tl2 



Ni(t) 



V d -l(a3 2 ) D 
V d (R 3 ) i 3 2 



N 3 (t) 



(C3) 



Vd(fl2) 



Vd-l(ai2) 1 Vd-l(a23) 



JV 2 (t). 



2. The triangle 

The triangular network studied here is depicted in 
Fig. EH (b) . All the containers C\ , C 2 and C3 are con- 
nected to each other and therefore the rate equations 
describing the time evolution of Ni(t), N 2 (t) and N 3 (t) 
have the same form. Equation (|C4|) is the rate equation 
that governs the dynamics in container C\. The corre- 
sponding dynamical equations for N%{£) and N 3 (t) are 
obtained in the same way but are not written down here 



C 2 and C3, that is Ni(t), N 2 (t) and N 3 (t) respectively, 
all have similar form, 

-w/o^^^^*') (C6) 

+A i4 {t-t')}, i = l,2,3. 

The middle container C4 has connections to all other 
containers C\ , C 2 and C3 . This leads to a rate equation 
for Ni(t) on the form 

Mt) = EJ=i^5^/o* / ^j(*')^4(*-*') 

x Ej=i V d {a A] )[K i3 {t - t') + A 4j (t - t')]. 

(C7) 

The rate equations derived in this subsection are used in 
the study of the junction, shown in Fig. ^] 



= Y V^fo dt '* f 2(t')*2l(t-t>) 



-%iter /o <ft'A/i(f JM* - «0 + a 12 (* - f )] 



- Y wm 1 Jo dt'W)[K 13 (t - f) + A 13 (t - f )]■ 

(C4) 

The behavior of Ni(t), N 2 (t) and A^(t) are very sen- 
sitive to changes in the network structure. The response 
from these changes are not easily predicted by Eq. I|C4() . 
Instead Eq. ijSlH) can be used to state a simplified ver- 
sion of Eq. I|C4|) . Note that this equation only is valid 
in the large time limit if the number of particles in the 
tubes is small. Applying Eq. I|59|l to container C\ in the 
triangular network gives 



JV x (t) = 



V d -i(a 2 i) D 



V d {R 2 ) 

D 

V d (Ri) 



N 2 {t) 



V d (R 3 ) t-31 6X ' 



Vd-i(ai2) , Vd-jXais) 



(C5) 



N!(t). 



Corresponding first order rate equations for N 2 (t) and 
N 3 (t) are found in the same way. 



3. The T-Network 

The T-shaped network under investigation in this sub- 
section is depicted in Fig. EH (c). Since all containers C\, 
C 2 and C3 are connected to container C4 and not to each 
other, the rate equations governing the dynamics in C\, 



APPENDIX D: ELIMINATION OF JUNCTIONS 
IN THE ASYMPTOTIC REGIME 

The network studied in this section is a generaliza- 
tion of the one depicted in Fig. 1151 Here, a junction 
having an arbitrary number of connections is studied 
and it will be demonstrated that the dynamical equa- 
tions governing the transport in a system having junc- 
tion points can be simplified in the large time limit. 
First, it will be shown that the current of particles in 
and out of a junction point eventually will balance each 
other out and that this occurs faster than the time it 
takes before an equilibrium particle distribution is at- 
tained throughout the network. This derives from the 
fact that the volume of the junction (proportional to a d ) 
is much smaller than the volumes of the containers, 
(proportional to Rf , i = 1, . . . , M), where d is the di- 
mensionality, a is the tube radius and M is the number 
of containers. The tubes connecting the junction are as- 
sumed to have the same radius. Second, it will be demon- 
strated how to simplify the transport equations in such 
a way that the dynamical variable for the junction N(t) 
can be eliminated completely. This might be desirable 
when working with large networks. 

Consider a junction with dynamical behavior con- 
tained in Nj(t) and volume Vd(a) that has connections to 
M containers with volumes Vd(Ri), i — 1, . . . , M. The 
equation governing the dynamics of the junction point 
written in Laplace space reads 



1 (a) 



V d (Ri) U 'J 



(«). 



(Dl) 
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This equation is a generalization of Eqs. i|29|) (see Sec- 
tion EJ- The junction point studied here is such that 
a <C Ri and since 



~v d Ji 



VARi) ~ Cd a { R, 



(D2) 



where q 



-r 



,/ i L (|) [r(^)] _1 7r _1/2 [found by using 
the formula for ad- dimensional sphere, see Section ITT) 
the term proportional (a/Ri) d can be neglected and a 
closed expression for Nj(s) can be obtained. Solving 
Eq. 1D1I) after eliminating the second term leads to 



«(i + ?E^iM») + AiiW] 



(D3) 



The inverse Laplace transform of Nj(s) is a sum of ex- 
ponentials 



P =i 



(D4) 



where a p is the residue of Nj(s) at pole s = —f3 p , 
(3 P G 5i > 0. The poles are given by the zeros of 



('a 



M 

E 



Djj cosh 



aiij sinh 



(D5) 



where qf^ 



slfjDi-j. When the currents in and out of 



the junction point balances each other out there is no 
accumulation of particles and Nj(t) = 0. This is easily 
verified from Eq. (|D4|I by evaluating the derivative in 
respect to t at t = oo 



can be used which was found for a two-node network. 
If the containers have volumes proportional to R d then 

^network ~JJ ("^") • If ^> d then Tj unc tion ^"network- 

This is is supported by the numerical calculation shown 
in Fig.lltjlwhcrc the curve corresponding to N^t) flattens 
out relatively fast. 

For networks where there are many junctions and con- 
tainers involved it might be desirable to reduce the num- 
ber of variables in the dynamical equations governing the 
particle transport. This can be done in the regime where 
Nj(t) = is valid and Eq. I|59|l is applicable. Applying 
Eq. H59fl for a junction point j having M connections and 
using Nj (t) — leads to 



M 



= Y,C l3 V d - l {t 



D, 



i=i 



V d (Ri) V d (a) 



(D7) 



In this way it is possible to express Nj(t) in terms of 
Ni(t), . . . ,NM(t). Inserting the solution to this matrix 
equation in the dynamical equations eliminates the ex- 
plicit dependence of Nj(t). Eqs. I|D8|) - (|D10|I shows this 
explicitly for the case of a three way junction depicted in 
Fig. 1151 where N4 (t) has been taken away completely 



I24D Vd-i(a) 
V d (R 2 ) 



N 2 (t) 



e 3i D v d 



(t 2i +l 3 j)D V d -i(a) AT / 1\ 
Z 5 V d (Ri) 



(D8) 



Mi) = 



V d (Ri) + V d (R 3 ) iV 3W 

P TW JV2(t) 



(D9) 



j t Nj{t) = Y,(-P P h P e- 0pt - as t 



(D6) 



p=i 



Let 



'junction 



be an estimate of the time it takes until 



Nj(t) = 0. It is related to the decay exponents (3 P (see 



Eqs. HD and which are zeros of Eq. i|D5|) . From 
Eq. 1D5I) it is clear that the zeros scales with S. (it is 
assumed that Dij ~ D and £^ ~ t) which leads to the 

Let rework be an estimate of 



estimate t; 



junction 



D ■ 



the time it takes to reach an overall equilibrium particle 
distribution in the network. A s a rough estimate Eq. 15' 



Mt) 



P V d (R 2 ) Nl "> + ~W V d (R 2 ) JS ^ t ) 

(eii+e 24 .)D v d _i(n) AT /,n 

P V d (R 3 ) iV 3W 

*14&4 + &4&4+4l4&4. 



(D10) 
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FIG. 1: Schematic picture of an arbitrary network built from 
containers and tubes. 



(a) 




(b) 



x = 

FIG. 2: Panel (a): A spherical compartment connected to a 
cylindrical infinitely long tube. If the tube radius is assumed 
to be small the compartment can be treated as ideally mixed 
at all times simplifying the dynamics in the container. The 
transport in the tube is reduced to a one dimensional diffusion 
problem. Panel (b) illustrates these simplifications. 
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FIG. 3: Behavior of f(y,z,t) defined in Eq. J7J for the two 
dimensional system shown in Fig. [panel (a), a/R = 0.1] at 
three different instants of time ti < t% <tz- The graph verifies 
the assumption that f(y,z,t) is approximately constant. In 
this two dimensional case f(y,t) — > l/2o as t — > oo. 
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FIG. 4: Schematic picture of a two-node network. 
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FIG. 5: Time dependence of the rate coefficients A(t), K(t) 
and a(t) from Eq. @. 
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FIG. 6: The solution of Eq. (I35H (solid line) is compared with 
a solution of Eq. I|55|l (dashed line). The solution of Eq. 145^ . 
given in 14911 . is represented by the dotted line. The curves 
decaying and growing describe Ni(t) and N2(t) respectively. 
The network structure is shown in inset. The volumes of the 
tube and the containers are equal and aft = 0.05. The initial 
distribution of particles is N 1 (0) /N tot = 1 and N2(0)/N to t = 
0, where N to t = iVi(0) + Ni(Q), indicated by the shading in 
inset. This figure clearly shows that Eq. 14911 does not lead to 
the correct values for Ni(t) and N2(t). The agreement with 
Eq. 155H is much better. 




FIG. 7: The natural logarithm of [Ni(t) - iVi(oo)]/iVtot for 
the three cases illustrated in Fig.HJ The labeling of the curves 
is the same as in Fig. E] The linear behavior is the evidence 
of the single exponential decay of the number of particles in 
the container C\. The slope gives the value of the decay 
exponent. The slopes of the solid and dashed lines are close 
to each other showing that Eq. I|55|l is capable of estimating 
the decay exponent well. The slope of the dotted line differs 
significantly from the others which illustrates that Eq. I49H 
can not describe the dynamics in an adequate way. Since the 
value of iVi(oo) is not the same as in all three cases [compare 
Eq. JSU and J^J] the value of Ni (0) -Ni (oo) will be different. 
This explains why all three curves do not coincide at t = 0. 
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FIG. 8: Dependence of the geometrical factor q 2 on the 
tube volume. (The volumes of the containers are equal 
Vd(Ri) = Vd(R2) = Vd-) q\ and r" 1 are related through 
r -1 = —q(D/£ 2 . The numerical solution of Eq. 1421 . which 
gives the exact values for q 2 , is represented by the solid line. 
It is compared with the values for q 2 ^ a (Eq. 1501 . dashed line, 
and q 2 b (Eq. 1561 . dotted line. The dotted line deviates sig- 
nificantly from the solid one as Vtube/Ki increases while the 
dashed line follows the exact solution better. This indicates 
that q 2 b provides a good estimate of the decay rate, even for 
large tube volumes. q\^ a can only be used for very small values 

Of Vtubc/Vd- 
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FIG. 9: A numerical verification of the assumption of well 
stirred containers. For simplicity reasons a cubic geometry, 
as depicted in panel (a), was chosen since the explicit shape of 
the container looses its importance for large reservoir volumes. 
The edge length is set to 2R where R is the radius of an 
equivalent spheric container. Panel (b) shows a numerical 
solution to the 2D-diffusion equation (solid line) compared 
to a numerical solution to Eqs. (dashed line). The cases 
are chosen so that a/R varies in three orders of magnitude 
showing increasing validity of the assumption of ideally mixed 
containers as a/R decreases. The initial distribution in all 
three cases are skewed (a delta function in the bottom right 
corner of the left container) which becomes important when a 
is large. For small a, r m i x is a lot shorter than Target and the 
skewed initial distribution will have time to smear out before 
particles start exiting and the shape of the initial distribution 
has no effect. 
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FIG. 11: The transport properties of the structure depicted in 
Fig. llOl fa'l. see Section lVlill fcase 1) for discussion. The curves 
depict a numerical solution of Eqs. ijCl^ and l|C2^ given in Ap- 
pendixEI] The curve depicting the number of particles in the 
middle container C 2 (solid line) has a maximum. This kind of 
behavior does not exist for the two-node network where there 
is only exponential growth and decay (see Fig. The net- 
work parameters were set to £12 = £23 = i, a/i = 1/5, a/Ri = 
1/10, a/Ri = 1/15, a/R 3 = 1/20. The initial distribution of 
particles is Ni(Q)/Nt„t = 1, N 2 (0)/N tot = N 3 (0)/N tot = 0, 
shown graphically in inset. 
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FIG. 12: The control of the particle arrival time for three 
different cases shown in inset: £ = £\ (solid line), £/£ 2 = 1/5 
(dash-dotted line) and ^£3 = 1/10 (dashed line). The system 
is otherwise equivalent to the one studied in Fig. 1111 The 
initial distribution of particles is shown graphically in inset. 



FIG. 10: networks used for case studies in Section IVHII Pan- 
els (a)-(d) corresponds to cases 1-4. 
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FIG. 13: The transport properties for the triangular network 
shown in Fig. 1101 (bt. see Section fyilll fcase 2) for discussion. 
The curves depict a numerical solution of Eq. I IC4I given in 
Appendix IH2l The curve depicting the number of particles 
in C 2 (solid line) has a minimum. This does not occur in the 
transport dynamics for the two-node network where there is 
only exponential growth and decay (see Fig. |SJ . The param- 
eters were set to t = i 2 i, t/a = 1, £ 23 / 'i = 1/10, £ 3 i/£ = 50, 
Ri/e = Ball = Ra/i = 4. The initial distribution of parti- 
cles is Ni(0)/N tot = 1, JV 2 (0)/JV to t = 0.5, N 3 (0)/N tot = as 
indicated by shading in inset. 
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FIG. 15: The transformation from a four-node network to a 
network involving a three way junction. 
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FIG. 14: The transport properties for the T-shaped network 
shown in Fig. 1101 (c), see Section I Villi (case 3) for discus- 
sion. The curves depict a numerical solution of Eqs. IC6I 
and {C7j given m Appendix |C2l The fi gure shows two ex- 
tremum points in the curve depicting the number of particles 
in container C 2 (solid line). It is not possible to have this 
kind of behavior for any of the cases studied involving two 
and three containers. The additional extremum point was 
produced by adding an additional container C4 to the middle 
container C 2 in the linear structure shown in Fie. 1101 fa). The 
network parameters used were £\ 2 = £, l/a = 3, £/£ 23 = 1/2, 
£/£ 2 4. = 10, Ri/i — 2. The initial distribution was set to 



Ni(0)/Ntat = 2/3, N 2 (0)/N tot = 0, AT 3 (0)/iV to( 
N4(0)/N to t = 1/3 according to shading in inset. 
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FIG. 16: The transport properties for the network involving 
a three way junction shown in Fig. 1101 (d), see Section IVIUI 
(case 4) for discussion. The curves depict a numerical solution 
of Eqs. and [with V d (R 4 ) = V d (a)\. The system pa- 
rameters were set to t\ 2 ja — 2, £13 /a = 10, £14 /a — 20, 
Ri/a — 2 and R 2 /a = R^/a — 4. Initial distribution 
of particles Ni(0)/N tot = 1, N 2 (0)/N tot = N 3 (0)/N tot = 
A r 4(0)/A r tot = is shown in inset. 
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FIG. 17: Structure of the network studied in Section I Villi 
case 5. The parameters describing the geometry are labeled 
in the same way as in Fig. 1101 The parameters used were 
a/l = 1/3, £ 23 /£ = 1, ks/i = 9, £ 43 /£ = 10, £ 35 /£ = 8, 
£ 36 /£ = 2, £ 28 /£ = 0.1, Ivxll = 1, £ 32 /£ = 1, Ri/a = 4, 
R 2 /a = 2, R 3 /a = 1, R 4 /a = 2.5, R 5 /a = 3, i? 6 /a = 2, 
#7 /a = 2, R s /a = 4. 
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FIG. 18: The solution of Eq. ^ for the network depicted 
in Fig. 1171 The panel includes graphs showing the transport 
dynamics for three different choices of initial distribution of 
particles. Panel (a): 7Vi(0) = N 2 (0) = N 3 {0) = N 4 (0) = 
AT 5 (0) = 0, jV 6 (0) = N 7 (0) = N 8 (0) = 1/3; Panel (b): 
N 3 (0) = N 5 (0) = JV 6 (0) = 0, N 2 (0) = N s (0) = 1/8 and 
ATi(O) = iV 4 (0) = N 7 (0) = 1/4; Panel (c): N 3 (0) = N 5 (0) = 
N 6 (0) = 0, JVi(0) = iV 7 (0) = 1/8 and N 2 (0) = iV 4 (0) = 
N 8 (0) = 1/4 where N,(t) = N t {t)/N tot , i = 1,...,8. The 
initial conditions are also illustrated in the insets: a darker 
shading indicates that more particles are injected into the 
container. 



